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ABSTRACT 

Aims. Deriving accurate frequencies, amplitudes, and mode lifetimes from stochastically driven pulsation is challenging, more so, 
if one demands that realistic error estimates be given for all model fitting parameters. As has been shown by other authors, the 
traditional method of fitting Lorentzian profiles to the power spectrum of time-resolved photometric or spectroscopic data via the 
Maximum Likelihood Estimation (MLE) procedure delivers good approximations for these quantities. We, however, show that a 
conservative Bayesian approach allows one to treat the detection of modes with minimal assumptions (i.e., about the existence and 
identity of the modes). 

Methods. We derive a conservative Bayesian treatment for the probability of Lorentzian profiles being present in a power spectrum 
and describe an efficient implementation that evaluates the probability density distribution of parameters by using a Markov-Chain 
Monte Carlo (MCMC) technique. 

Results. Potentially superior to "best-fit" procedure like MLE, which only provides formal uncertainties, our method samples and 
approximates the actual probability distributions for all parameters involved. Moreover, it avoids shortcomings that make the MLE 
treatment susceptible to the built-in assumptions of a model that is fitted to the data. This is especially relevant when analyzing solar- 
type pulsation in stars other than the Sun where the observations are of lower quality and can be over-interpreted. As an example, we 
apply our technique to CoRoT* observations of the solar-type pulsator HD 49933. 
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1. Introduction 

Our understanding of the Sun's structure has been revolutionized 
over the last three decades by the application of helioseismology, 
where the surface manifestation of acoustic modes (p-modes) 
are used to probe the interior. Seismology is now being applied 
to stars using precise rapid photometry from space and high- 
precision radial velocity measurements from the ground. Low 
degree p modes were identified in several main-sequence and 
sub-giant stars (see e.g., Bedding & Kjeldsen 2006). 

Due to the lower signal to noise ratio of the data obtained 
from stars, and the more poorly constrained stellar models, there 
is a greater risk of getting things wrong. It is critical to a mean- 
ingful analysis to have a measure of the reliability and accuracy 
of the observed stellar frequency spectrum. Even a few incor- 
rectly identified frequencies, i.e., frequencies that are not intrin- 
sic to the star but an artifact of the data processing or of instru- 
mental origin, will significantly complicate the analysis and may 
even lead to a misinterpretation of the observations. Indeed, it is 
a case where quality overrides quantity, that is, a few well ob- 
served modes frequencies are more useful in model fitting than 
a large number of poorly determined or unreliable frequencies. 

Here we describe a Bayesian approach to extracting frequen- 
cies from a noisy stellar oscillation spectrum. The Bayesian 



* The CoRoT (COnvection, ROtation and planetary Transits space 
mission, launched on 2006 December 27, was developed and is operated 
by the CNES, with participation of the Science Programs of ESA, ESAs 
RSSD, Austria, Belgium, Brazil, Germany and Spain. 



method treats the problem of identifying frequency peaks in 
terms of probabilities. The probability that a peak in the spec- 
trum is an oscillation mode is calculated using basic physical 
assumptions about the shape of the modes profile and the na- 
ture of the background noise. If the data warrants it, additional 
prior knowledge can be included to increase the complexity of 
the mode shape, aiding in mode identification and the analysis 
of rotational effects (i.e., how frequency is split by rotation). The 
actual application of the Bayesian approach to compute the pos- 
terior probability distributions of models evaluated using simu- 
lated or real data is computationally complex and can quickly 
exceed our computational resources. Therefore, we use an im- 
plementation of the Markov-Chain Monte Carlo (MCMC) tech- 
nique, which significantly reduces these demands. The applica- 
tion of Bayesian treatments to solar-type oscillations is not new 
(e.g., Brewer et al. 2007). The implementation of Markov-Chain 
Monte Carlo techniques for a Bayesian analysis of Lorentzian 
profiles has also already been described (Benomar 2008). 

Because we are very concerned about the reliability of the 
data (i.e., whether individual peaks in the spectrum are intrin- 
sic or spurious) and specifically want to avoid interpreting noise, 
however, we advocate a more conservative approach. In particu- 
lar, we use the mode height parameter of our models so that we 
can quantify the probability that a frequency peak is real or sim- 
ply due to noise. This eventually allows us to arrive at a model 
that is just complex enough to comply with the data and its noise 
properties. 
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In the next section we describe the application of Bayes the- 
orem to the problem of identifying stellar modes in an oscillation 
spectrum. In section 3 we describe our application of the MCMC 
technique to the computations. In section 4 we test our method 
on simulated data and in section 5 we apply our methodology 
to the recently obtained observations from CoRoT for the star 
HD49933. 



2. The probability of a p-mode power spectrum 
model 

In the following subsections we introduce the Bayesian for- 
malisms including the likelihood function and the definition and 
role of priors. But first we review the basic properties of solar- 
type oscillations (e.g. (see e.g., Appourchaux et al. 1998). 

From the equation of a damped harmonic oscillator forced by 
a random function, the average power of the Fourier transform 
of the displacement can be approximated as, 



(P(co)) 
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where Pj is the power of the Fourier transform of the forcing 
function, and ojq and y are the angular mode frequency and the 
damping term, respectively. If the spectrum of the forcing func- 
tion varies slowly with the frequency, the resulting average p- 
mode profile can be approximated by a Lorentzian profile as 



P(f,v,h,r) 
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where v is the mode frequency, / is some frequency value in the 
spectrum, and h and r\ are the height and width (77(7-) = (mY x , 
with t being the mode lifetime) of the profile in the power den- 
sity spectrum, respectively. For non-radial modes (t > 0) stellar 
rotation will split the mode into a 2( + 1 multiplet. This can be 
modeled by adding 2( Lorentzian profiles, spaced with integer 
multiples of the rotation frequency around the central mode. In 
this case, Equ. 2 is replaced by 



P(f,v,v tot ,h,T)= J] — 
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where v ml is the rotational splitting. The mode heights h m of 
the additional profiles with m + are determined by the ob- 
served geometry, specifically, the inclination angle. In general, 
the mode height is related to the mode amplitude according to 
a 2 = Jirjh. Usually a Maximum Likelihood Estimation (MLE) 
algorithm is used to find a best fit between the observed power 
spectrum and a model like 
P m (f) = b + Z i P(f,v i ,h i ,T i ), 

with b being the background noise (which is assumed to be lo- 
cally white) and P( f, v,, hi, r ( ) being the power of the individual 
profiles at some frequency (-bin) 

2.1. Application of Bayes Theorem 

The Bayes' Theorem can be stated as follows: 



p(A\D,D = 



p(A\I)p(D\AJ) 
P(D\I) ' 



(4) 



where p(A\D, I) is the probability for some hypothesis A given 
the data D and the prior information /. Note that / is also called 
the posterior probability of A. p(A\I) is the prior probability of A, 



p(D\A, I) is the likelihood function, and p(D\I) is the global like- 
lihood. The latter serves as a normalizing constant and ensures 
the total probability summed over all hypotheses equals one. An 
extensive introduction to this field and the relevant methods can 
be found in Gregory (2005). 

In order to apply the Bayesian formulation to any problem, 
one has to identify and assign probabilities for the individual 
terms in Equ. 4. In the following subsections we discuss our 
choices. 

2.2. The Likelihood Function 

The statistic of a power spectrum is given by a x 2 distribution 
with 2 degrees of freedom. The probability density that a model 
spectrum m matches an observed power spectrum o corresponds 
to 



p(m) = Y\ 
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where P (f) is the observed power at frequency /, and P m (f) de- 
notes the corresponding expectation value, i.e., the power given 
by the model (e.g., Toutain & Appourchaux 1994). Equ. 5 is 
used in the MLE method as a likelihood function and has al- 
ready been suggested for a Bayesian treatment of solar-type os- 
cillations (Appourchaux 2008; Benomar 2008). The Bayesian 
formalism, though, can contain additional prior probabilities for 
each of the parameters from which the model spectra are con- 
structed. 

2.3. The role of the priors - improving the MLE approach 

The simplest model of a p-mode profile consists of four param- 
eters: mode frequency, mode height, mode lifetime, and back- 
ground offset. Physically, the mode height and mode lifetime 
parameters are correlated. However, here we refrain from using 
this information and assume them to be independent. This en- 
ables us to use the mode height parameter as a device to distin- 
guish between a detection and a non-detection of p-mode signal. 
We define the following two priors, called ignorance priors be- 
cause they make no assumptions about the physical properties of 
the object being modeled: 

1 . For a mode with frequency v, mode lifetime t, background 
offset b we assume a uniform prior of the form 



P(x\l) = 



1 



(6) 



where x stands for the respective parameters v, r, and b. That 
is we assume that the probability is equal across the range of 
parameters and makes obvious the nomenclature, ignorance 
prior. For the prior of the frequency parameter this is an 
appropriate choice, since we cannot a priori exclude, for an 
observed spectrum, any value within the range defined by 
the upper and lower p-mode frequencies. The mode lifetime 
and the background offset, as long as they do not vary across 
the fitted part of the spectrum by more than a magnitude, 
are mostly determined by the global properties of the power 
spectrum and therefore also warrant a uniform prior. 

2. For a mode with height h we adopt a modified Jeffrey's prior 
of the form 



P (h\D = 



1 
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where 



2.4. Treatment of rotational effects using ignorance priors 



h = log 



t h' min + h max \ 

V ^m'm ' 



(8) 



and h' ■ (> 0) expresses the "strength" of the prior. Lower 
values of h' min decrease the probability associated with h max . 
In contrast to the uniform prior, the Jeffrey's prior is nor- 
mally employed to ensure a uniform probability density per 
decade range (e.g., the prior probability between 0.001 and 
0.01 is the same as between 0.01 and 0.1). The modified prior 
behaves like a Jeffrey's prior for h > h' min , yet allows for val- 
ues smaller than h' min . In this regime it acts like a standard 
uniform prior and avoids a logarithm of 0. 

With the Bayesian framework restricted to ignorance priors, any 
evidence for significant power at some frequency / can be han- 
dled by this combination of priors. The modified Jeffrey's prior 
of the mode height provides a solution for the problem that in 
real observations power due to noise can be mistaken for ac- 
tual intrinsic signal. But by allowing the likelihood function to 
deliver significantly higher probabilities for mode heights well 
above the noise level, it reduces the mode height prior's pref- 
erence for a non-detection. The choice of h! in this context 

mm 

is therefore nothing more than the expression of an "odds ratio 
condition". Consider the posterior probability of the mode height 
at some value h, without considering the additional parameters, 
which is 



p(h\D,I) oc p(h\I)p(D\h,I). 
The ratio of the probabilities 
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is an odds ratio for the models "mode height = h" and "mode 
height = h' min ". It can also be seen as a Bayesian weighting of 
a likelihood ratio, or as a strength-of-evidence indicator similar 
to a SNR. O p r i or is the prior odds ratio, while OnkeUhood is of- 
ten called the "Bayes factor". The following equation is then an 
"odds ratio condition"/or the definite detection of the Lorentzian 
profile: 



O 



1 
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o„ 



<o 



likelihood ■> 



(ii) 



If O cond is larger than On ke uhood, the model "mode height = h' min " 
is favoured, and there is not enough evidence for the detection of 
a Lorentzian profile with a height greater than this threshold. For 
example, O CO nd = 10 5 requires a likelihood function value of a 
given h that is at least 10 5 times that of h' min in order for the pro- 
file to considered as detected. That this condition arises for the 
mode height parameter is quite intuitive, since the mode height is 
the parameter that determines whether a Lorentzian profile rises 
above the noise level. In order to ensure a constant odds ratio 
condition for a data set, independent of the number of modes to 
be detected, the geometric mean of the mode height priors can 
be used when more than one Lorentzian is evaluated. 

To summarize, the mode height clearly is a scale parameter, 
rather than a location parameter (see Gregory 2005). Choosing 
h' ml „ allows one to set a lower limit for mode peak detection. The 

mm ' 

MLE procedure does not implicitly allow for such a restriction 
and, thus, makes it prone to overfitting. 



The effects of stellar rotation complicate the picture tremen- 
dously. If rotational splitting is suspected, one can replace sin- 
gle Lorentzian profiles with a sum of profiles. The central 
Lorentzian of the rotationally split mode can be treated just like 
a single profile, i.e., it has the same parameters and correspond- 
ing priors. The number of suspected split components, however, 
depends on the spherical degree of the mode. The components 
themselves are characterized by the rotational splitting and their 
mode heights, which relative to each other and the central peak 
depend on the stars inclination angle (see, e.g., Gizon & Solanki 
2003) If the SNR conditions are very poor and/or a prelimi- 
nary mode identification is not possible, we chose not to use 
the inclination angle as a parameter but to set the heights of 
the split components as individual free parameters (see Equ. 3). 
Each Lorentzian profile in the model can a priori be "equipped" 
with a number of rotationally split components conforming to 
the highest value of ( expected to be present in the data. The 
corresponding height priors, then, only allow significant mode 
heights of those components in the rotationally split multiplets 
for which there is enough evidence in the data. This way, no 
preliminary mode identification is necessary. Alas, using this ap- 
proach the number of parameters needed for the model increases 
tremendously. For high SNR data, where the degree of the modes 
becomes even visually apparent through the rotational splitting, 
a preliminary mode identification is certainly a more sensible 
approach. In such a case, though, our method is not needed any- 
way. 



3. Application using Markov-Chain Monte Carlo 

The analytic evaluation of complex models in a large parame- 
ter space in terms of Bayesian probability soon reaches the limit 
of computer resources. Stochastic methods can provide a suf- 
ficient sampling of the parameter regions of interest at a much 
smaller cost. In particular applications of the MCMC-technique 
have therefore gained increasing popularity and the (Bayesian) 
problems to which MCMC nowadays are applied range from the 
detection of planets (Gregory 2007), to the analysis of solar-type 
pulsators (Brewer et al. 2007), to spot modelling of active stellar 
atmospheres (Croll 2006). 

Basically, MCMC performs a sequential walk through the 
parameter space of a specific problem. The MCMC technique 
probes solutions to the Bayesian equations by generating a ran- 
dom set of model parameters at some point and then progressing 
through parameter space via a biased random walk. This bias di- 
recting the chain of steps through parameter space is provided 
by the Bayesian probabilities themselves. Specifically, each pa- 
rameter of the model is incremented or decremented by some 
random fraction of a predefined step width and the procedure 
then accepts or declines this combination of steps. The condition 
of acceptance is provided by the Metropolis-Hastings algorithm 
(Hastings 1970), which is based on the ratio of probabilities of 
subsequent parameter configurations before and after the step. 
In order to comply with this algorithm, the random steps do not 
necessarily have to be sampled from a symmetric proposal distri- 
bution. However, using such a distribution (uniform, Gaussian, 
etc.) is more intuitive and facilitates the necessary calibration 
of the algorithm (see below). MCMC scales with an increasing 
number of parameters and it is an ideal tool for solving our prob- 
lem. We note that the Metropolis-Hastings algorithm evaluates 
only the relative probability between different parameter config- 
urations. The normalisation factor of Equ. 4, which is often very 
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difficult to evaluate, does not need to be calculated and will be 
dropped here after. 

3. 1. Semi-automated MCMC calibration 

In order for the MCMC-implementation to operate efficiently 
the acceptance rate must be properly calibrated. When there are 
more than two parameters it can be shown that the acceptance 
rate should be about 0.25 in order to minimize correlations be- 
tween the different parameters (Roberts et al. 1997). Because the 
acceptance rate depends mainly on the step width of each pa- 
rameter, the step widths, themselves, have to be calibrated. The 
calibration process becomes more and more cumbersome as the 
number of parameters increases. We, therefore, implemented the 
following automated MCMC calibration scheme, similar to that 
described in Gregory (2005). 

Given an observed power spectrum, several non-overlapping 
frequency windows are defined that envelope the range of the in- 
dividual mode frequency parameters. The number of Lorentzian 
profiles to be investigated per window is specified, and the lower 
and upper constraints for the mode height(s), the mode life- 
time^), etc. are defined. The model spectrum is initialized with 
appropriate starting parameters, such as equidistant frequencies 
within the frequency windows and mean values for the mode 
height, mode lifetime, and so on. The MCMC algorithm is 
started with a step width of 



f (-^) — 0.1(x mflJC %>min) 



(12) 



for each model parameter x, where x max and x min are the up- 
per and lower limits already used to define the prior probabil- 
ities. During this "burn-in" phase our MCMC algorithm eval- 
uates, for every single step and for every single parameter, the 
relative probability of the models according to Equ. 4, using 
the described likelihood function and priors. It will take several 
hundreds to thousands of iterations to approach the parameter 
regions of maximum likelihood. When the "burn-in" phase is 
completed the model should be close to the global maximum of 
probabilities. For the next few thousand iterations, the individ- 
ual acceptance rate for each parameter is evaluated every hun- 
dred or so steps. Simultaneously, the MCMC step width of the 
respective parameter is slowly adjusted to approach a desired ac- 
ceptance rate, which depends on the total number of parameters 
in the model. For example, we found that in order to achieve a 
combined acceptance rate of about 0.25 for a model consisting 
of ~ 70 parameters, the individual acceptance rates had to be set 
to roughly 0.94. 

Once the desired individual acceptance rates remain fairly 
stable, our algorithm switches to the standard MCMC routine, 
where the probability of a new configuration, and therefore the 
probability of its acceptance, is evaluated after all parameters 
have been slightly changed according to a random walk. At this 
point the MCMC routine is set to automatically test the parame- 
ter space. 

3.2. Evaluating the MCMC results 

Since the acceptance probability for each step is calculated us- 
ing the Bayesian posterior probability, the MCMC routine will 
compute distributions (marginal distributions) for each involved 
parameter. After a sufficient number of iterations the marginal 
distributions for all model frequencies, mode heights, mode life- 
times, and background offsets can be estimated from histograms 
that count how frequently respective values of each parameter 



were visited. An estimate for the validity of the tested model 
can be obtained by examining these distributions. If any strong 
asymmetries, multiple local maxima, and other obviously non- 
Gaussian features appear they can be used to guide one toward 
an improved model. Regardless, the uncertainties for all param- 
eters can be obtained from the cumulative distribution function 
of the normalized histograms. 



4. Simulations 
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Fig. 1. Results of the Bayesian analysis for one of the 100 simulated 
data sets. Upper panel: (a) the power density spectrum calculated from 
the simulated data (dashed line) and the most probable model spectrum 
derived from the analysis (solid line); (b) the marginal distributions 
of the mode frequency parameters for both Lorentzian profiles of the 
model. Middle panel: the corresponding marginal distributions of the 
mode height parameters. Dashed and solid lines refer to the respective 
distributions presented in the top panel (b). Bottom panel: the marginal 
distribution of the mode lifetime parameter. 



In order to test our method with parameters similar to a typi- 
cal CoRoT run we generated 100 time series of 60 days duration 
each following the procedure in Chaplin et al. (1997). Each sim- 
ulated data run represents a different realisation of the signal de- 
scribed in Tab. 1 . We calculated the power spectra of each time 
series and applied our Bayesian analysis. The frequency range 
was between 130 and 140 d _1 . Two Lorentzian profiles were fit- 
ted, each with adjustable mode height but equal mode lifetime. 
The upper and lower limits for all parameters, which also enter 
the Bayesian formalism via the prior probabilities, are presented 
in Tab. 2. h' . was set to 10~ 6 , and h max to 0.15, and the geomet- 
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Table 1. Input parameters for the simulated solar-type pulsator. 





/i 


h 


frequency 


132.3 


136.6 


rms amplitude 


1 


1 


estimated mode height 


0.033 


0.033 


mode lifetime 


0.5 d 


0.5 d 



white noise model with SNR of 0. 1 in the time domain 



Table 2. Parameter limits for the analysis of the simulated solar-t 
pulsator. According to Equ. 7, h' min was set to 10~ 6 . 





min 


max 


y[d-'] 


130 


140 


h 





0.15 


T[d->] 


0.1 


1.5 


b 


0.001 


0.01 



ric mean of the mode height priors was used. This correspond 
an odds ratio condition of ~ 10 5 for the maximum value h max ; 
to a condition of ~ 10 4 for the actual input mode heights. r 
results are based on 300000 iterations of the MCMC routine. 



4.1. Evaluation of the simulations 

Fig. 1 shows a typical power spectrum calculated from one 
the 100 test data sets, along with the most probable fit to me 
data. All input parameter values fall within the borders of the 
derived marginal distributions. We present the median values for 
all the parameters (but see discussion in Gregory (2005)) and, 
importantly, the uncertainties. A closer inspection of Fig. 1 re- 
veals that not all the derived parameter values fall within the 
lcr uncertainties of the parameter values of the simulated data. 
Indeed, all the parameters for this example are underestimated: 
The median of the f\ distribution deviates from the input value 
by 2.99cr. The corresponding mode height is underestimated by 
0.37cr. f 2 is underestimated by 0.63<x, with the corresponding 
mode height deviating by 1 .2<x. The mode lifetime also deviates 
by 0.84<x. Such deviations are expected and are a direct conse- 
quence of the stochastic nature of the signal. 

Fig. 2 shows the results for all 100 simulations. The me- 
dian values derived from our analysis are compared to the sim- 
ulation input values and are plotted against the corresponding 
normalized probabilities of these input values, evaluated via the 
marginal distributions. As expected, the input values are not al- 
ways reproduced to within lcr (p > 0.682), but are consistently 
found within 3<x. The scatter of the derived median values of 
the parameters around the input values behaves as expected. The 
median of the frequency parameter is shown to be distributed 
symmetrically, while the mode lifetime and mode heights follow 
a log-normal distribution. The figure also illustrates the overall 
accuracy with which each parameter can be reproduced. The fre- 
quency parameter shows a lower accuracy limit of about 0.1%. 
The corresponding lower limits for the mode height and mode 
lifetime parameters are much larger at about 40%. 

In addition, Fig. 2 shows cumulative histograms for the pa- 
rameter value probabilities of the evaluated simulations. These 
histograms indicate that the scatter of the probabilities of the 
input values, for individual realisations, can be roughly approx- 
imated by a normal distribution. As the derived marginal distri- 
butions for all parameters provide a consistent picture, we argue 
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Fig. 2. Right panels: distribution of the derived median mode parameter 
values (v, h, r) compared to the input values. Positive/negative values 
indicate under-/overestimation of the parameters. The open squares and 
plus signs in the top and middle panels identify the results for the two 
Lorentzian profiles in the model. Left panels: cumulative histograms for 
the probabilities of the simulation input values determined from all 100 
simulations. 



that the probability densities obtained from a single data set can 
be trusted. 

4.2. The importance of the mode height prior 

As was argued in Sec. 2.3, the mode height prior in its given form 
allows one to adjust the fitting to better account for noise. The 
simulations described in the previous Section provide us with an 
excellent example for this rationale. In Fig. 3 we show one of 
the simulations that includes by chance an additional power ex- 
cess in the frequency range of consideration due to noise and/or 
the stochastic nature of the input signal. Using these data we 
tested how looking for of a (non-existent) third Lorentzian pro- 
file would be influenced by the h' mjn parameter of the mode 
height prior. We evaluated the probabilities for the existence of 
3 Lorentzian profiles in the spectrum and performed the analy- 
sis with two different values of h' . = 10~ 8 and h! . = 0.01. 

mm mm 

The weaker prior with h' ■ = 0.01 failed to correctly identify 
the noise peak as noise. This is similar to the behavior of the 
MLE method, or any Bayesian application that does not im- 
plicitly treat the problem of over-fitting due to, e.g., noise. The 
strong prior with h' mjn = 10~ 8 , correctly distinguished between 
the noise peak and the two real profiles (see Fig. 3) and, in fact, 
determined a maximum in probability for a mode height of zero 
for the third. 

In the case of the strong prior, the marginal distribution of 
the frequency parameters for profile 2 and profile 3 are good ap- 
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Fig. 3. 3 Lorentzian profiles are fitted to simulated data that contain 2 frequencies. Top panel: power density spectrum of one of the 100 simulations 
mentioned in Sec. 4. The 2 input frequencies are marked by thick dashed lines. Additional power at ~ 131.7d~' is either due to noise or, more likely, 
due to the stochastic nature of the signal. Middle panels: marginal distributions of the frequency parameters calculated after 300000 iterations. 
Each panel contains the results for one of the three fitted profiles. The thick lines show the results for h' mjn = \Q i , the thin dashed lines indicate 
h' min = 0.01. Lower panels: corresponding marginal distributions of the mode height parameters for each of the Lorentzian profiles and for both 
values of h'. (10~ 8 and 0.01). 



proximations to a normal distribution and match the two input 
frequencies. Indeed, profile 2 is not influenced by the additional 
power excess in its vicinity. Although the marginal distribution 
shows a maximum probability for the frequency parameter at the 
location of the additional power excess, the marginal distribution 
also assigns significant probability densities to the whole param- 
eter range. But more to the point the marginal distributions for 
the height prior, which again are well-shaped for profile 2 and 
profile 3, show a most probable mode height of zero for profile 
1. That is, the profile 1 is unnecessary to fitting this range of 
frequencies with only two profiles are clearly detected. 

In the case of the weak prior, the marginal distributions of 
the frequency and mode height parameters of profile 1 and pro- 
file 2 intersect and mix. Hence, there seems to be evidence for 
a third Lorentzian profile which is comparably strong to the ev- 
idence for the actual Lorentzian profile at 132.3 d _1 . Moreover, 
profile 1 and profile 2 cannot be easily separated. Thus, their 
mode parameters cannot be unambiguously determined. In sum- 
mary, with a weak mode-height prior (mimicking MLE meth- 
ods) a third mode is detected where none exist and the near by 
mode information is distorted. 



5. Application to the CoRoT target HD 49933 

The question of detecting Lorentzian profiles is only relevant 
for data sets, where the SNR is high enough to see evidence 
for solar-type pulsation, but low enough to be questionably 
tractable with methods used in helioseismology. The ground- 



braking CoRoT data have all the qualities necessary for the tech- 
nique to work: 

1. good data quality devoid of strong instrumental signal or 
aliasing 

2. clear indication of Lorentzian profiles at least in the region 
of maximum pulsation power 

3. increasingly poor S/R further away from the region of maxi- 
mum pulsation power. 

With the application to such data, we will demonstrate the prac- 
tical utility of our Bayesian method. We applied our technique 
to the 60 days of CoRoT N2-data of HD 49933 obtained dur- 
ing the Initial Run from February to March, 2007. A detailed 
description of the data extraction and reduction can be found in 
the CoRoT book (see Baglin et al. 2006, and references therein). 
The data set consists of more than 163 000 measurements sam- 
pled each 32 s. A detailed summary of the stellar properties and 
the CoRoT observations is given by Appourchaux et al. (2008). 

Intrinsic stellar noise-like signal, mainly due to turbulent 
convection, generates significant power in the frequency range 
of pulsation and effects the accurate determination of mode pa- 
rameters. For the Sun Harvey (1985) modeled the background 
signal with a sum of power laws. Aigrain et al. (2004), and 
more recently Michel et al. (2008), use P(v) = 2;Pj> with 
Pi = fl,^r,/(l + (2wt,-v) c '), or hereafter P, = A,/(l + (fi,v) c <), 
to model the solar granulation signal, where v is the frequency, 
t, is the characteristic times scale of the ith component and C, is 
the slope of the power law. The normalisation factor a, is cho- 
sen such that = J Pj(v)dv corresponds to the variance of the 
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Fig. 4. Top panel: power density spectrum of HD 49933. Light grey: 
original power density spectrum; dark grey: the original spectrum 
smoothed with a 20 //Hz boxcar average; thick solid line: global fit ac- 
cording to Equ. 13; dotted line: global fit without the Gaussian term; 
dashed lines; white noise and the 3 power law components of the global 
fit. Middle panel: enlargement of the top panel. Dotted line: background 
fit plus the Lorentzian profiles fitted to the residual power density spec- 
trum. Bottom panel: original power density spectrum with the granula- 
tion signal removed. Solid line: Gaussian term of the global fit. 



stochastic signal in the time series. Whereas the slope of the 
power laws was arbitrary fixed to 2 in Harvey's original model, 
Appourchaux et al. (2002) were the first to suggest a different 
value. Recently, Aigrain et al. (2004) and Michel et al. (2008) 
have shown that, at least for the Sun, the slope is closer to 4. 
Each power law represents a different physical process with time 
scales for the Sun ranging from minutes in the case of granula- 
tion to days in the case of stellar activity. 

Appourchaux et al. (2008) used a single power law for the 
background signal in the frequency range of pulsation and fitted 
the corresponding parameters simultaneously with the /?-mode 
parameters to the observed power density spectrum. We have 
separated the analysis of the background signal from the analy- 
sis of the pulsation signal, in which case one has to model the 
pulsation signal in the power density spectrum with a dedicated 
term in the global background fit. Kallinger et al. (2008) have 
shown that the power excess hump due to pulsation can be ap- 
proximated by a Gaussian function. Hence, our global model to 
fit the heavily smoothed power density spectrum of HD 49933 
consists of a superposition of white noise, three power law corn- 



Table 3. Global fit parameters. The amplitudes of the power laws (A,), 
the height of the Gaussian part (P g ), and the white noise components 
(P„) are given in ppm 2 //iHz. The power law time scales (B,) are given 
in seconds and the center (v max ) and width (<x) of the Gaussian part are 
in /jHz. One-sigma error estimates are given in brackets. 





Ai 


A 2 


A 3 




Bi 


B 2 


B 3 


Power laws 


25743 (6200) 


22 (6) 


2.2 (0.11) 




252311 (25000) 


16477 (1950) 


1637 (62) 






^max 


cr 


Gaussian term 


0.500 (0.015) 


1657 (28) 


538 (26) 


White noise 


P„ 


= 0.33 (0.005) 





ponents, and a power excess hump approximated by a Gaussian 
function: 



P(v) = P, 



(=i 



(B,v) 4 



-(v-v„,„) 2 /(2(r 2 ) 



(13) 



where P„ is the white noise component, A, and B, are the ampli- 
tudes and characteristic time scales of the power laws, P g , v max , 
and cr are the height, central frequency, and width, respectively, 
of the Gaussian term. The resulting global fit (and its compo- 
nents) is shown in the top and middle panel of Fig. 4 along with 
the original and heavily smoothed power density spectrum. The 
corresponding parameters of the global fit are given in Tab. 3. 
The fit without the Gaussian component (dotted line in top and 
middle panel of Fig. 4) enables us to estimate the local noise and 
is used to separate the background signal form the pulsation sig- 
nal. The residual power density spectrum rescaled to the white 
noise level, shown in the bottom panel of Fig. 4, is used for the 
subsequent analysis. 

5.1. MCMC analysis 

The frequency range between the lower and upper limits be- 
yond which we see the pulsation signal disappearing in the noise 
was subdivided into 16 windows. To each, a chosen number of 
Lorentzian profiles was fitted, based on our Bayesian algorithm. 
The model parameters are presented in Tab. 4. This subdivision 
into windows was only chosen to accelerate the burn-in phase. It 
has no influence on the results, as long as the resulting marginal 
distributions for the mode frequency parameters are not inter- 
sected by the window borders. Yet, it still enables us to perform 
a global fit which includes the influence of Lorentzian profiles 
also from adjacent windows. 



Table 4. Parameter limits for the analysis of HD49933. The height prior 
parameter was set to h' min = 10~ 8 . 





min 


max 


v[pHz] 


1219.9 


2618.1 


h [ppm 2 /pHz] 





6 


r[d->] 





1.5 


b [ppm 2 /pHz] 


0.001 


0.5 



For the first analysis we decided to consider only 2 
Lorentzian profiles in each window. Moreover, we also kept the 
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mode lifetime uniform for all Lorentzian profiles, expecting that 
the marginal distribution of this parameter would tell us whether 
this assumption is valid. 

Fig. 5 shows the power density spectrum, on which our anal- 
ysis is based, together with the fitting windows we defined, 
and the most probable model derived after ~ 1.3 million iter- 
ations. The model manages to reproduce the observations quite 
well, and there remain only very few frequency ranges where 
a visual inspection seems to indicate additional power. Beyond 
v > 2450 /Aiz several Lorentzian profiles are fitted with fairly 
well defined frequencies, but with a probability maximum for a 
mode height of zero. Two profiles, listed in Tab. 5 as P25 and 
P26, indicated in the figure by dashed lines, have ambiguous 
frequency distributions ranging over the whole fitting window. 
What is shown is a biased choice to fit the general picture. 
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Fig. 6. Top panel: marginal distribution of the mode frequency param- 
eter of one of the profiles fitted to the CoRoT data. The median value, 
and the lower and upper ltr-values are indicated by the dashed lines. 
Middle panel: same as top panel but for the corresponding mode height 
parameter. Bottom panel: same as top panel but for the mode lifetime 
parameter for a simultaneous fit to all Lorentzian profiles. 
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Fig. 7. Top panel: (a) the fitted region in the power density spectrum, 
where the p modes have the highest SNR. The most probable model 
(thick solid line) is presented in comparison to the observations (dot- 
ted line). The borders of the fitting windows are indicated by the thick 
dashed lines, (b) The marginal mode frequency distributions for 3 pro- 
files per window. Six profiles are identified (black), while the additional 
profiles are ambiguous across the whole frequency range (gray lines). 
Frequencies identified by Appourchaux et al. (2008) are shown as thin 
dashed lines. Bottom panel: The corresponding mode height distribu- 
tions. 



5.2. On the evidence fort = 2 modes 

We repeated our analysis for the region with the highest SNR, 
but included for each fitting window a third Lorentzian pro- 
file in our model. The results are shown in presented in Fig. 7. 
The probability distributions for the additional profiles consis- 
tently favours a mode height of zero. Moreover, the correspond- 
ing mode frequency distributions only vaguely contain regions 
of higher probability. We therefore conclude that, based on our 
conservative approach, any additional power in the power den- 
sity spectrum is due to noise or to the stochastic nature of the 
clearly detected modes. From our perspective, focussing on the 
"detection" of profiles, the data do not present convincing evi- 
dence for I = 2 modes. 



Fig. 6 presents an example for the marginal frequency distri- 
bution and mode height for one of the Lorentzian profiles. As ex- 
pected, the former can be very well approximated by a Gaussian, 
while the latter follows a log-normal distribution. The bottom 
panel shows a quite narrow and smooth mode height distribution 
for a simultaneous fit to all Lorentzian profiles, indicating a con- 
stant mode lifetime. The data therefore do not seem to warrant 
different mode lifetimes. It is thus save to assume that any varia- 
tions among the various profiles for this parameter are accounted 
for within the uncertainty given by the mode lifetime distribu- 
tion. Surprisingly, we do not find well-defined multi-modal dis- 
tributions for any of the frequency parameters involved. Some 
ambiguities in the distributions are only apparent at very high 
frequencies and with the poorest SNR. Additional profiles, ei- 
ther due to modes of higher degree or rotational splitting, still 
could be present, if the frequencies are very well separated. 



5.3. On the evidence for rotational splitting 

The failed detection of additional profiles, and the single-mode 
nature of the frequency distributions of detected profiles, sug- 
gests that the effects expected from rotational splitting can 
most likely be neglected in the analysis of the CoRoT data of 
HD 49933. Nonetheless, in order to perform a more rigourous 
test, we analysed the same frequency region shown in Fig. 7 in- 
cluding these effects. Our model used 2 Lorentzian profiles per 
fitting window, each of which contains 2 additional features cor- 
responding to rotational splitting of ( = 1 modes (see Equ. 3). 
Accordingly, a new global parameter, the rotation frequency, was 
introduced. Appourchaux et al. (2008) report the detection of 
a particular spectral feature at 3.4 //Hz which apparently corre- 
sponds to the HD49933's rotation frequency. We therefore used 
this value as a reference and allowed a scatter of about 0.2 //Hz 
in either direction. This range corresponds to twice the classical 
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Fig. 5. Observed power density spectrum (gray) and fitted Lorentzian profiles. Vertical lines indicate the borders of the fitting windows (see text). 
The dashed lines in the bottom panel show two profiles for which our Bayesian technique results in ambiguous frequency distributions ranging 
over the whole fitting window 



frequency resolution (Ar) _1 of the CoRoT time series, where 
AT is the length of the data set, and was used to define the bor- 
ders of a uniform prior. The heights of the rotation profiles were 
implemented as usual, using modified Jeffrey's priors for these 
parameters. We expected a well defined distribution for the ro- 
tation frequency to emerge from the analysis. However, we find 
no preferred value for this parameter, as is shown in Fig. 8. The 
distribution is mostly consistent with sampling noise. 

Concerning the mode heights, alternating detections and 
non-detections of the rotation profiles would indicate a dif- 
ference in their central profiles between radial and non-radial 
modes. We cannot find any such differences, since the mode 
height distributions for the rotation profiles are consistent with a 
null result. While the star certainly does rotate, it seems that the 
signal is too close to the noise level to support any need for these 
parameters. In other words, rotational effects, if present, seem to 
be ill-defined and do not need to be considered for the frequency 
analysis of this data set. As an example of a near ideal case, the 
bottom panel of Fig. 8 shows corresponding results from a sin- 
gle mode using 60 days of data of the Sun observed as a star by 
the VIRGO (green band) instrument on board the SOHO space- 
craft (see Frohlich et al. 1997, and references therein). Fitting 
only the small range between 3090 and 3110 yuHz and using 
the same profile model as for HD 49933, the rotational splitting 
becomes apparent. 




3.4 3.45 

v ra >Hz] 




0.4 0.45 

v ra >Hz] 



Fig. 8. Top panel: the marginal distribution of the rotational splitting 
parameter derived as explained in Sec. 5.3. The thick black line is a 
boxcar average using 10 data points. Bottom panel: the same but for the 
analysis of VIRGO data of the Sun. 
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Fig. 9. Echelle diagram showing the identified frequencies (solid error 
bars). Values belonging to profile 25 and 26 (see Tab. 5) are indicated by 
the dashed error bars. Frequencies found by Appourchaux et al. (2008) 
are displayed in grey. 



5.4. Results 

As a consequence of the results from the previous two sections, 
we use only 2 modes per window without any rotational split- 
ting for our global fit. For 6 of the 32 fitted profiles, the iden- 
tification is ambiguous (or the most probable mode heights are 
close to zero) due to the low SNR. All other profiles are clearly 
detected and their parameters have smooth, single-mode distri- 
butions. The results of our analysis are presented in Tab. 5. The 
lcr-errors have been derived from the cumulative distribution 
functions, which are evaluated using the marginal distributions 
for all parameters. We find that the data are consistent with a 
mode lifetime of roughly 0.5 days which does not appear to vary 
significantly across the spectrum. The confidence for this over- 
all result can be calculated via the mode height parameter. As we 
have used h' min = 10~ 8 , the result gives an odds ratio condition 
O CO nd = 4.07 x 10 6 . 

Therefore, the obtained parameter distributions are at least 4.07 x 
10 6 times more probable than a model spectrum that only con- 
sists of the background offset. The profiles 25 and 26 are only 
listed for sake of completeness. Due to the ambiguity in their 
marginal frequency distribution, apparent in their given uncer- 
tainties, we do not assign credibility to these values. The final 
Echelle diagram is displayed in Fig. 9. The frequencies estimated 
as I = 1 by Appourchaux et al. (2008) agree very well with our 
corresponding values, and the lcr-uncertainties are comparable 
in both studies. The remaining frequencies, however, differ due 
to the cited authors' explicit assumption of I — 2 modes, which 
we do not detect. 

Although we deliberately did not include fixed mode height 
ratios for modes of different degree in our model, the resulting 
mode height ratios are mostly consistent within the uncertainties 
with a fixed ratio. The obvious outlier in Fig. 10 at ~ 2280yuHz 
is due to profiles 25 and 26, which we already marked as highly 
suspicious. 



Table 5. Results of the CoRoT data for HD 49933. Profiles P 25 and P 
26 are ambiguous detections. Profiles P 29 to P 32, which are not listed, 
have a most probable mode height of zero. 
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6. Discussion 

With high quality data and evidence for additional and/or closely 
spaced modes, such as might appear for nonradial modes split 
by rotation, we believe that it might be advantageous to change 
the peak Lorentzian model by combining the corresponding 
Lorentzian profiles and fitting them as a group. Eventually, once 
the quality of the data becomes even better, we think that a 
more classical Bayesian treatment replacing the ignorance priors 
with specific information from theory or previous observations, 
is more appropriate because the potential gain in information 
outweighs the danger of over-fitting. Benomar (2008) has per- 
formed such an analysis for the data set also presented in this 
paper, and obviously obtained different results. It thus remains 
questionable whether the Initial Run CoRoT data of HD 49933 
is already good enough for such an approach. 

6.1. On the global likelihood of p-mode profiles 

In this paper we have used Bayes' theorem to solve a parameter 
estimation problem. One of the biggest advantages of Bayesian 
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Fig. 10. The mode height ratios for subsequent modes listed in Tab. 5 
(e.g., hpi/hpi) are shown as filled circles. The dotted lines correspond 
to the 68.2% confidence limits. The dashed line indicates a height ratio 
of unity. The outlier at ~ 2280 yuHz is due to profiles 25 and 26, which 
are ambiguous detections. 



analysis, however, is to perform model comparison. This is 
achieved by first calculating the global likelihood via integra- 
tion over all model parameters. The global likelihoods are sub- 
sequently used to compare the different models themselves. It 
would be interesting to use parallel tempering (e.g., Benomar 
2008), an algorithm which allows the Markov chain to more eas- 
ily sample the complete parameter space, to perform such model 
selection by calculating the global likelihood of a model via inte- 
gration over the tempering parameter (see, e.g., Gregory 2007). 
However, we think that the validity of applying this kind of 
model selection to p-mode analysis needs to be carefully tested 
first. 

As a first trial run, we have investigated how the global like- 
lihood of three simple models behaves in a region of pure noise. 
We studied the region between 4600 and 4700 ^Hz in the power 
spectrum of HD 49933 (see Fig. 1 1) which is far beyond the re- 
gion of pulsation. At such high frequencies, the power spectrum 
should be dominated by the noise properties of the instrument. 
Following the description given in Gregory (2005) on how to 
obtain the global likelihood, we evaluated three models 

- Ml: constant noise level, 

- M2: constant noise level + 1 Lorentzian profile, and 

- M3: constant noise level + 2 Lorentzian profiles 

The noise level was allowed to vary between 0.001 and 0.5 
ppm 2 /iHz _1 . All Lorentzian profiles were allowed a mode life- 
time between and 1.5 days, mode heights between and 6 
ppm 2 /iHz _1 , and frequencies between 4600 and 4700 pHz. We 
used 10 chains with decreasing but equidistant tempering param- 
eter. In order to study the effects of the likelihood function on the 
gobal likelihood we neglected all priors in the calculations. The 
results are best expressed using the odds ratios 



4600 



4620 



4640 4660 
frequency [uHz] 



4680 



4700 



p(M\\D,I) 
p(M2\D, I) 



exp(-240) and 



p(M2\D, I) 
P(M3\D,I) 
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Depending on the choice of the tempering parameters used in 
the parallel tempering process, as well as on differences in the 
numerical integration scheme required to calculate these num- 
bers, the end result varies, but the overall magnitudes of dif- 
ference in probability remain the same. The likelihood function 
alone seems to prefer models with (more) Lorentzian profiles, 
even if there are none (related to solar-type oscillations) present 
in the data. 



Fig. 11. The high-frequency region devoid of p-mode signal used for 
testing the global likelihood comparison. Note that the scale is similar 
to Fig. 5. 

This is an issue demanding clarification. It will be required to 
test the influence of the allowed parameter ranges in order to see 
where this preference arises in the integration over the param- 
eter space. Moreover, the influence from different kind of prior 
probabilities on the global likelihoods also needs to be studied. 
To conclude, our initial test suggests that we cannot yet trust the 
global likelihood completely when comparing models of power 
spectra (to real data at least) until we are able to clarify this issue. 

7. Conclusion 

We have presented a conservative approach (using only igno- 
rance priors) for a Bayesian analysis of solar-type p modes with 
a minimum of parameters. The method uses a mode-height prior 
that allows us to more easily and reliably distinguish real peaks 
from those that arise stochastically from noise, in reference to a 
threshold set by the user. We show how the marginal distribu- 
tions of all the parameters can be obtained. Our tests with sim- 
ulated data succeeded by reproducing the simulation input val- 
ues and not being confused by randomly occurring noise peaks. 
Furthermore, we have applied our method to the Initial Run 
CoRoT data of HD 49933, where we can easily and unambigu- 
ously identify at least 26 p modes. We fail to detect 1-2 modes 
or effects produced by rotation. We also approached the possi- 
bility of performing Bayesian model comparison in the analysis 
of solar-type p modes. First results suggest, however, that more 
tests need to be performed in order to understand how model 
comparison works in a power spectrum. Finally, we would like 
to point out a companion paper (Kallinger et al. 2009). Therein, 
we show that frequencies extracted in this work are in excellent 
agreement with model frequencies of a solar-calibrated model 
that coincides with the spectroscopically determined position of 
HD 49933 in the H-R diagram. 
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